﻿# -*- coding: utf-8 -*-
# Code sous licence CeCill
# http://www.cecill.info/licences.en.html
# This is a license compatible with the GNU GPL license
# and that conforms to the principles of French law.
# Authors: X. FABREGES, Nicolas Bruyant
"""
TDO datas treatment
===============
Provides
    1.Read files from gage acquisition card
    2.Perform fast frequency demodulation
    3.Optionnal adjacent averaging
    4.Integrate pick up signal


TODO:
-----
    Split in function
    Create unit tests
    

"""
__version__ = "Version 0.1, SVN-$Revision: 14 $"
"""
Revision allow a good tracking of data processing  errors
Refers to SVN revisions: DO NOT edit manually
"""


from scipy import int16, mean, arange, zeros, real, fft, argmax, pi, savetxt, fromfile, r_, ones, convolve, exp, floor, cos
from scipy.fftpack import fftfreq
import sys,os
from pylab import *


#from  ..devlib.importfilters import *
from  importfilters import *

#import cProfile

# Declaration globale des slices (faster)
slice1 = [slice(None)]
slice2 = [slice(None)]
slice1[-1] = slice(1, None)
slice2[-1] = slice(None, -1)


def smooth(x, window_len=11):
    """
    Perform adjacent avaraging

    Parameters
    ==========
        x : datas to filter
        window_len : averaging length
    Returns
    =======
        y Filtered datas

    """
    s=r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]

    w=ones(window_len, 'd')
    y=convolve(w / w.sum(), s, mode='valid')
    return y


def local_trapz(y):
    # Methode des trapezes modifiee : constante supprimee
    """
    Trapezoidal integration optimized for speed

    Parameters
    ==========
        y: datas to integrate

    Returns
    =======
        integrated datas

    Remark
    ======
        Use slice1 et slice2 global variables to improve speed

    """
    global slice1
    global slice2

    return ((y[slice1]+y[slice2])).sum(axis=-1)


def tdo_fft(inputfile, outputfile, WindowType='Cos', time_le = 4096):
    """
    Perform Fourier transform and return frequency evolution

    Parameters
    ==========
        inputfile, outputfile : filenames
        WindowType: 'Rect', 'Cos', 'Bess'
        time_le:  timewindow length

    Returns
    =======
        integrated datas

    Remark
    ======
        Use slice1 et slice2 global variables to improve speed

    TODO
    ====
        Split in Two functions: file import and data treatment

    """

    # Input parameters
    fourier_le = 1024   # Fourier length
    dfmin      = 1      # Frequency resolution
    dt         = 2e-8   # timestep of acquisition = 0.6sec/3e7 points

    # Lecture du fichier
#    fid = open(inputfile, 'rb')
#    fid.seek(512)   # Skip useless header
#    V = fromfile(fid, int16, -1, '')
#    fid.close()
    V = loadfileGage(inputfile)

    pstart     = 1      # First timewindow
    pend       = int(floor((len(V)-fourier_le)/time_le))  # Last timewindow

    V = V-mean(V)  # Remove zero frequency contribution to FT
    t = arange(0, fourier_le)*dt

    # Approximation of main frequency
    Vf         = abs(real(fft(V[0:fourier_le])))
    tf         = fftfreq(fourier_le, dt)
    fmax       = zeros((pend+1-pstart, 6))
    fmax[0, 1] = tf[argmax(Vf[0:int(fourier_le/2)])]
    fmax[0, 0] = time_le*dt/2

    #load PU datas
    filepath, filename = os.path.split(inputfile)
    dbdt_data = loadfileGage(os.path.join(filepath, 'p' + filename[1:]))
    #calculate field
    field = field_integration(dbdt_data,
                             method='end',
                             offset_calculation_length=2E4)
    # adjust field to remove drift
    #field=field_adjust(field,method='end')
                         
    # keep only mean field value for each slice
    SlicedField = zeros((pend-pstart+1))
    for i in xrange(pstart, pend):
        fmax[i, 3] = np.mean(field[i*time_le:i*time_le+fourier_le])
        fmax[i, 4] = np.min(field[i*time_le:i*time_le+fourier_le])
        fmax[i, 5] = np.max(field[i*time_le:i*time_le+fourier_le])
    
    

    # Calculation of constants
    expon   = -2j*pi*t
    deltaf0 =  tf[1]/1000


    # Windowing
    if  WindowType=='Cos':       
        # Cosinus
        window = arange(0,fourier_le)
        window = 1-cos(window*2*pi/(fourier_le-1))
    elif  wWindowType=='Bess':    
        # Kaiser-Bessel
        window  =  arange(0,fourier_le)
        window  =  0.402-0.498*cos(window*2*pi/(fourier_le-1))+0.098*cos(window*4*pi/(fourier_le-1))+0.001*cos(window*6*pi/(fourier_le-1))
    else: # WindowType='Rect'
        # Rectangular
        window  =  ones(fourier_le)

    # Precise determination of oscillating frequency via DFT
    for i in xrange(pstart, pend):
        print "Calculation for slice ", i
        # Utilisation de la dernière valeur comme point de depart
        a = fmax[i-1, 1]
        V_temp = window*V[i*time_le:i*time_le+fourier_le]

        # Previous frequency spectral weight
        # Complex exponential time consuming ! 
        #TODO Need a smarter way to perform this calculations
        deltaf = deltaf0
        essaimax=np.abs(local_trapz(V_temp*exp(expon*a)))

        # Calculation of local derivative of Fourier transform
        # If derivative positive, then search for frequency in growing direction
        if abs(local_trapz(V_temp*exp(expon*(a+deltaf)))) > essaimax:
            while abs(deltaf) > dfmin:
                F = abs(local_trapz(V_temp*exp(expon*(a+deltaf))))
                if F > essaimax:
                    essaimax = F
                    a += deltaf
                else:
                    deltaf = -deltaf/5

            # Store frequency
            fmax[i, 0:2] = [(i*time_le+fourier_le/2)*dt, a-2.5*deltaf]
        # Lower frequency otherwise
        else:
            while abs(deltaf) > dfmin:
                F = abs(local_trapz(V_temp*exp(expon*(a-deltaf))))
                if F > essaimax:
                    essaimax = F
                    a -= deltaf
                else:
                    deltaf = -deltaf/5

            # Store frequency
            fmax[i, 0:2] = [(i*time_le + fourier_le / 2)*dt, a + 2.5 * deltaf]


    # Save calculation of smoothed datas in file
    fmax[:, 2] = smooth(fmax[:,1])[5:-5]

    savetxt(outputfile, fmax)
    return fmax

i=0
while i<len(sys.argv):
    # Scans to extract
    if sys.argv[i]=='-i':
        inputfile=sys.argv[i+1]
    elif sys.argv[i]=='-o':
        outputfile=sys.argv[i+1]
    i += 1



res = tdo_fft(inputfile,outputfile,time_le = 1024)
#plot frequencyvs field
plot(res[:,3],res[:,2])
plot(res[:,4],res[:,2])
plot(res[:,5],res[:,2])